Computationally Efficient Direction-of-Arrival Estimation Algorithms for a Cubic Coprime Array

In this paper, we investigate the problem of direction-of-arrival (DOA) estimation for massive multi-input multi-output (MIMO) radar, and propose a total array-based multiple signals classification (TA-MUSIC) algorithm for two-dimensional direction-of-arrival (DOA) estimation with a coprime cubic array (CCA). Unlike the conventional multiple signal classification (MUSIC) algorithm, the TA-MUSIC algorithm employs not only the auto-covariance matrix but also the mutual covariance matrix by stacking the received signals of two sub cubic arrays so that full degrees of freedom (DOFs) can be utilized. We verified that the phase ambiguity problem can be eliminated by employing the coprime property. Moreover, to achieve lower complexity, we explored the estimation of signal parameters via the rotational invariance technique (ESPRIT)-based multiple signal classification (E-MUSIC) algorithm, which uses a successive scheme to be computationally efficient. The Cramer–Rao bound (CRB) was taken as a theoretical benchmark for the lower boundary of the unbiased estimate. Finally, numerical simulations were conducted in order to demonstrate the effectiveness and superiority of the proposed algorithms.


Introduction
The multiple-input multiple-output (MIMO) radar [1][2][3] has received much attention in array signal processing, due to its capability to enhance spatial resolution, improve parameter identifiability, and achieve more degrees of freedom (DOFs) than traditional phased array radars [4]. Recently, a new MIMO system, called massive MIMO, has attracted great interest both in industry [5] and academia [6][7][8]. Compared to the traditional MIMO geometry, massive MIMO systems are equipped with a large number of antennae at the base station (BS), which brings a number of advantages: for example, higher throughput, improved spectral output, and enhanced link reliability [9][10][11][12][13].
Base stations rely on massive MIMO systems to uplink sounding signals, and compute channel knowledge to enable MIMO beam-forming [14][15][16]. Therefore, direction-of-arrival (DOA) estimation is important for massive MIMO], as it is needed to complete the predicted capacity gains [17,18]. Various DOA estimation algorithms have been proposed in recent years, such as improved multiple signal classification (MUSIC) [19], estimation of signal parameters via rotational invariance techniques (ESPRIT) [20], and the propagator method (PM) [21]. However, all these classic algorithms are conducted with a uniform array [22][23][24][25]. In a uniform array, inter-element spacing is limited to half a wavelength to avoid the problem of phase ambiguity [26].
In recent years, sparse arrays, such as nested [27] and coprime [28] arrays, have received considerable attention due to their capacity to improve estimation accuracy, enhance DOF, and mitigate the mutual coupling (MC) effect [29][30][31]. Practically speaking, a nested array is composed of at least one uniform linear array. For a uniform array, because the inter-element spacing is no larger than half a wavelength to avoid phase ambiguity, severe MC effects still occur. To reduce the MC effects, a general coprime linear array (CLA) which contains two uniform linear subarrays was proposed by [28]. The number of sensors of these two subarrays are defined as M and N. The inter-element spacing of the sensors is set to Nλ/2 and Mλ/2, respectively, where λ represents the wavelength of the carrier wave. The authors of [28] verified that a CLA with M + N − 1 sensors can achieve O(MN) DOFs and improve the performance of the DOA estimation, where O(·) means the DOFs. A spatial smoothing technique-based algorithm was proposed in [31] which constructs a positive semi-definite covariance matrix to apply the subspace-based method. However, this algorithm involves a large number of snapshots to vectorize. In [32], a decomposed algorithm was proposed for CLA which separately processes the two subarrays by applying the conventional MUSIC algorithm and conducting a global angular spectral search. Meanwhile, the phase ambiguity problem is solved by combining the two estimates of two subarrays. A search-free algorithm was devised to estimate DOAs for coprime arrays with low computational complexity [33].
In practical applications, two-dimensional (2D) DOA estimation contains more important information. Consequently, a large number of investigations have been conducted into 2D DOA estimation generation [34][35][36]. The authors of [34] presented a general coprime planar array, proposing a partial spectral search algorithm to estimate 2D DOA which uses a linear relation to decrease the complexity. The authors of [35] proposed a combined ESPRIT algorithm which can achieve auto-paired angle estimation with lower complexity. The authors of [36] employed a reduced dimensional method which involved transforming the 2D spectral peak search to one-dimensional (1D) one within a small search sector. The authors of [37] proposed an initialization-based parallel factor method which utilized the PM to obtain initial estimates, then used these estimates to reconstruct the signal matrix; in this way, the complexity was significantly decreased. However, the above-mentioned studies, [34][35][36][37], are based on spectral peak searching, which results in high computational complexity.
For a two-dimensional array design, such methods can obtain superior performance in terms of azimuth angle estimation, but the elevation angle estimation is less good due to the fact that the array elements are distributed over a plane. From the perspective of improving performance in terms of the 2D DOA estimation of elevation angle, we propose considering the three-dimensional array model. Furthermore, the traditional antenna array is limited by the half-wavelength limit on inter-element spacing that avoids phase ambiguity. As a result, the array aperture is limited, which leads to poor resolution. Especially in the massive MIMO environment and in the context of 5G communication, which requires more antennae, the traditional array structure-with its limited array aperture, poorer resolution, lower estimation accuracy, and more severe MC effect-is not viable for practical applications. As a result, it is of great significance to expand the one-dimensional linear array and two-dimensional planar array models and introduce a three-dimensional coprime cubic array structure. Such designs can be used in massive MIMO, 5G communication, and UAVs. A coprime cubic array (CCA) geometry-a configuration which has been proven to outperform conventional uniform cubic arrays in terms of DOA estimation [38,39]-was proposed in [40]. In this paper, we exploit the same type of CCA configuration.
For CCA geometry, existing approaches, such as subspace-based algorithms [34][35][36], process the subarrays separately. This approach only employs the auto-covariance matrix of the subarray, resulting in achievable DOFs being damaged and making it impossible to solve the ambiguity problem caused by inter-element spacing larger than half a wavelength. In this paper, we propose a total array-based MUSIC algorithm (TA-MUSIC) for 2D DOA estimation, where not only the auto-covariance matrix, but also the mutual covariance matrix, is employed. Additionally, based on [34], we verified that, by using the coprime property, ambiguous estimates can be eliminated. To alleviate the computational burden, we present an ESPRIT-based MUSIC (E-MUSIC) algorithm which employs the computationally efficient ESPRIT algorithm to obtain initial DOA estimates, then conducts 1D Sensors 2022, 22, 136 3 of 17 spectral peak searches. Finally, numerical simulation results are presented to demonstrate the effectiveness of the proposed algorithms with CCA geometry.
The main contributions of our research are summarized as follows: (a) We propose a TA-MUSIC algorithm with CCA geometry to enable massive MIMOs to fully employ DOFs and achieve superior DOA estimation by employing both the auto-covariance matrix and the mutual covariance matrix of the entire array. In addition, we verify that by using the coprime property, the proposed algorithm can suppress the ambiguity problem. (b) We propose an E-MUSIC algorithm for 2D DOA estimation, which can effectively decrease the complexity of the classic MUSIC algorithm. After utilizing the ESPRIT algorithm to initialize and obtain a rough estimation, we then conduct a fine search within a smaller sector to achieve lower complexity. (c) Our numerical simulation results confirm that the proposed algorithms outperform the classical ESPRIT algorithm and the PM algorithm in DOA estimation.
The remaining sections of this paper are organized as follows. The array geometry and signal model are described in Section 2. The proposed algorithms are presented in Section 3. The complexity analysis is presented in Section 4, and the advantages of the proposed system discussed. Numerical simulations are provided in Section 5, and Section 6 discusses the conclusionsr.
Notations: We use lower-case bold characters for vectors and upper-case bold characters for matrices. The characters (·) T and (·) H denote the transpose and the conjugate transpose, respectively. The characters and ⊗ represent the Khatri-Rao product and Kronecker product, respectively. The element diag(·) is a diagonal matrix which uses the elements of the matrix as its diagonal elements. The element E(·) is statistical expectation; D m (·) is a diagonal matrix in which the m-th row of the matrix is employed; and angle(·) and arctan(·) are the phase operator and the arctangent function, respectively.

Signal Model
In this paper, we employ a CCA configuration which can further enlarge the array aperture and improve DOA estimation performance [40].
As described in [40], a CCA configuration incorporates two uniform cubic subarrays. One subarray is set with M 1 × T 1 × P 1 sensors. The other is set with M 2 × T 2 × P 2 sensors. For subarray 1, the inter-element spacing is determined as d x1 = M 2 λ/2 for the x-axis direction, d y1 =T 2 λ/2 for the y-axis direction, and d z1 = P 2 λ/2 for the z-axis direction, where λ represents the wavelength. In subarray 2, the inter-element spacing is d x2 = M 1 λ/2, d y2 = T 1 λ/2 and d z2 = P 1 λ/2, for the x-, y-, and z-axes, respectively. The total number of the sensors is therefore computed according to the equation Figure 1 depicts an example of a CCA configuration where M 1 = T 1 = P 1 = 4 and M 2 = T 2 = P 2 = 3.
If we denote the outputs of two subarrays as per [34], then: where S = [s 1 , s 2 , . . . , s K ] T ∈ C K×L represents the signal matrix and L denotes the number of snapshots. N i ∈ C M i T i P i ×L (i = 1, 2) is the additive Gaussian white noise matrix. The variance and mean are defined as σ 2 n and zero, respectively.  , where k  represents the elevation angle notes the azimuth angle of the k-th signal;

2) sensors and is presented by
. Then, we take subarray i as an example to illustrate the pro rithm in the following part. Accordingly, the subarray has If we denote the outputs of two subarrays as per [34], then: represents the signal matrix and L denotes of snapshots.
( 1,2) is the additive Gaussian white noise matri ance and mean are defined as 2 n  and zero, respectively.
i  sensors and is pr The directional matrix can also be denoted as: Here, we use A x_i , A y_i , A z_i as directional matrices of subarray i. A x_i represents the x-axis, A y_i denotes the y-axis, and A z_i is the z-axis. The character denotes the Khatri-Rao product, and D m (·) denotes a diagonal matrix which employs the m-th row values of the matrix.

Proposed Method for DOA Estimation
In this section, we examine the classic MUSIC algorithm [32] for CCA geometry. We then propose a TA-MUSIC algorithm to solve the problems encountered with the classic MUSIC algorithm. Note that the TA-MUSIC algorithm involves spectral peak searching within the global angular sector, generating a high level of computational complexity; we therefore propose using an E-MUSIC algorithm to alleviate the computational burden. To be more explicit, we take subarray 1 as an example to explain the proposed algorithm.

Review
In [32], the authors ran the subarrays of CCA separately. Correspondingly, the covariance matrix for output signal of the subarray i can be acquired using the formula: Therefore, according to eigenvalue decomposition (EVD), we have: whereD s,i represents the largest K eigenvalues in the signal subspace andÊ s,i represents the corresponding signal subspace matrix. The noise subspace, which is composed of remaining eigenvalues, is denoted byD n,i and the corresponding noise subspace matrix isÊ n,i .
According to the 1D MUSIC spectrum [19], we can denote the 3D MUSIC spatial spectral function for subarray i as

TA-MUSIC Algorithm
Unlike the algorithms that conduct the subarrays separately [32][33][34], in the proposed TA-MUSIC algorithm, both covariance matrices of the total array are employed to estimate DOAs by stacking the outputs of the two subarrays, as follows: where T . From this, the covariance matrix of the total array can be estimated using the formula: whereÊ s is the signal subspace matrix andÊ n denotes the noise subspace matrix.D s and D n denote matrixes which contain the largest K eigenvalues and the remaining eigenvalues, respectively. Then, according to the orthogonality that exists between the signal subspace and noise subspace, we can denote the spatial spectral function as: , the TA-MUSIC spectrum can be obtained without any ambiguous phase problems arising. This is demonstrated in Lemma 1, below [34]: We suppose (θ k , φ k ) to be the real DOA for the k-th signal. In addition, only one DOA estimate pair(θ k ,φ k ) generates a spectral peak in the TA-MUSIC spectrum. Therefore,(θ k ,φ k ) is the estimates of the real DOA(θ k , φ k ).
As described above, the proposed algorithm enables spectral searching within the global angular sector, which results in expensive computational complexity. In the following section, we propose using an E-MUSIC algorithm to significantly reduce the computational complexity.

The E-MUSIC Algorithm
From Equation (1), we know that the subarray with M i × T i × P i sensors has uniform inter-element spacing and the properties of a Vandermonde matrix, indicating that the ESPRIT algorithm, which has high computational efficiency, can be employed to achieve initial DOA estimates.
Next, we rearrange the rows ofÊ s,1 to obtain a new signal subspace,Ê s,1 Similarly toÊ s,1 , we then segmentÊ s,1 intoÊ y1,1 andÊ y2,1 , denotingÊ s,1 (1 : (T 1 − 1) M 1 P 1 , :) andÊ s,1 (M 1 P 1 : T 1 M 1 P 1 , :), respectively. Then we perform EVD on Ψ y1 =Ê ÷ y1,1Êy1,2 to obtain Φ y,1 andv a k . Likewise, we can achieve Φ z,1 andŵ a k by applying the same algorithm. In this part of the calculation, in order to decrease the complexity, we chose to explore only two variants. In the process of computing Φ x,1 ,Φ y,1 , and Φ z,1 , estimated DOAs (û a k ,v a k ) with automatic pairing can be obtained using the same matrix, T [28]. According to Equation (19), if the inter-element spacing is larger than half a wavelength, the phase ambiguity problem arises and all the DOA estimates we obtain will be ambiguous. However, according to the coprime property, the phase ambiguity problem can be eliminated and unambiguous values can be attained, as demonstrated in Lemma 1. To obtain more accurate DOA estimates, the unambiguous values (û ini k ,v ini k ) are employed to initialize the spatial spectral function in Equation (13).
Unlike the TA-MUSIC algorithm, which generates high computational complexity due to its use of 2D spectral peak searching within a global angular sector of θ k ∈ (0 • , 90 • ), φ k ∈ (0 • , 180 • ), the proposed E-MUSIC algorithm transforms the 2D spectral peak search into a 1D spectral peak search over the sector, u k ∈ (−1, 1), v k ∈ (−1, 1), which is significantly more computationally efficient.
For the k-th signal, the initial estimate (û ini k ,v ini k ) is initially employed. From this, a new 1D spatial spectral function can be obtained: A search is conducted over a small sector w ∈ (0, 1), and a 1D spatial spectral peak search is then performed to achieve an improved DOA estimate for w k . Similarly, by using w k andv ini k , we can construct the 1D spatial spectral function of u as: From this, it is possible to obtain more accurate estimates from a 1D spectral peak search. The final estimates can be obtained by means of:

Detailed Steps
The detailed steps of the E-MUSIC algorithm are provided as follows.
Step 1. Compute the covariance matrixR of the total array according to Equation (8), then decomposeR to attain the covariance matrixR of subarray 1.

Computational Complexity
In this subsection, we compare the complexity of the proposed TA-MUSIC algorithm, the E-MUSIC algorithm, and the classic MUSIC algorithm [32] using the same CCA for each. For the TA-MUSIC algorithm, the total complexity is where L represents the number of snapshots, G = G 1 + G 2 , G 1 =M 1 T 1 P 1 , G 2 =M 2 T 2 P 2 , and τ denotes the searching step. For the E-MUSIC algorithm, the total complexity is For clarity, we list the computational complexity of the mentioned algorithms in Table 1, where M 1 = 4, T 1 = 4, P 1 = 4, M 2 = 3, T 2 = 3, P 2 = 3, L = 200, K = 2 and the searching step τ is set to 0.02. We can clearly see that the proposed TA-MUSIC algorithm has the almost the same complexity as the TSS algorithm, while the proposed E-MUSIC algorithm has the lowest computational complexity compared to the others.  Figure 2 compares the complexity with increasing number of sensors of the three algorithms, classic MUSIC, TA-MUSIC and E-MUSIC. From this chart, we can observe that the computational complexity of these three algorithms increases with the number of sensors. It can also be seen that E-MUSIC has the lowest complexity among the three. Figure 3 compares the complexity of the three algorithms as the number of snapshots increases. This chart indicates E-MUSIC's superiority due to lower complexity. Figure 4 compares the complexity under the variant searching step. It clearly indicates that the computational complexity increases as searching decreases.
sensors. It can also be seen that E-MUSIC has the lowest complexity amo Figure 3 compares the complexity of the three algorithms as the number of creases. This chart indicates E-MUSIC's superiority due to lower comple compares the complexity under the variant searching step. It clearly indi computational complexity increases as searching decreases.

Degree of DOF
The proposed TA-MUSIC algorithm can obtain 1  The algorithms in [34][35][36] divide the array into two subarrays, then esti separately. In this way, much of the array information is not fully employed, r a loss of DOFs. The method described in [34][35][36] can achieve a greater numbe as

Degree of DOF
The proposed TA-MUSIC algorithm can obtain M 1 T 1 P 1 + M 2 T 2 P 2 − 1 DOFs, which indicates TA-MUSIC can detect M 1 T 1 P 1 + M 2 T 2 P 2 − 2 target signals by using The algorithms in [34][35][36] divide the array into two subarrays, then estimate DOA separately. In this way, much of the array information is not fully employed, resulting in a loss of DOFs. The method described in [34][35][36] can achieve a greater number of DOFs, as DOF methodsin where Specifically, the E-MUSIC algorithm has M 2 T 2 P 2 DOFs, Thus, the proposed algorithms can detect many more signals than conventional algorithms with the same sensors. From Equations (31)-(33), we can describe the relationship of DOFs of the three algorithms as:

Advantages
The performance of the TA-MUSIC and E-MUSIC algorithms are summarized as follows:

•
The proposed TA-MUSIC performs better DOA estimations by employing all of the array information, including the auto-covariance matrix and the mutual covariance matrix, whereas E-MUSIC only utilizes the auto-covariance matrix information. In addition, TA-MUSIC can fully achieve DOFs of M 1 T 1 P 1 + M 2 T 2 P 2 − 1, while the algorithms in [34][35][36] only achieve M 2 T 2 P 2 DOFs. Furthermore, E-MUSIC can obtain M 1 T 1 P 1 DOFs, which are larger than M 2 T 2 P 2 (M 1 > M 2 ).

•
The proposed E-MUSIC algorithm significantly reduces computational complexity by transforming the 2D spatial spectral peak search into a 1D one, whereas the algorithms in [34][35][36] and TA-MUSIC require 2D spectral peak searching.

•
The proposed TA-MUSIC algorithm can attain paired angles automatically and outperforms the conventional ESPRIT and PM algorithms in DOA estimation performance.

Cramer-Rao Bound
In this part, we provide the Cramer-Rao bound (CRB) [41] derivation of DOA estimation with a CCA. The signal covariance matrix is achieved through R s = SS H /L.

Comparison of the DOA Estimation Performance of Different Algorithms
In this subsection, we present the RMSE performance of various al CCA, specifically, the PM [42] algorithm, the ESPRIT [43] algorithm, an TA-MUSIC and E-MUSIC algorithms. It is assumed that there are two directions of 1 1 ( , ) (10 ,20 )      and 2 2 ( , )    (20 ,30 )   . Figure 7 dep against the SNR (signal-to-noise ratio), where the number of snapshots i 8 depicts the RMSE against the number of snapshots, where SNR = 0dB. In be seen that estimation performance improved as with increasing SNR, a posed TA-MUSIC algorithm outperforms the others, because it utilizes information, from both the auto-covariance matrix and the mutual-covar addition, Figures 7 and 8 show that the E-MUSIC algorithm, which u scheme to estimate the signals' direction of arrival, can achieve promine tion performance, along with the high computational efficiency demonst 2 and 3. The root mean square error (RMSE) was exploited as the estimation performance metric, defined by the following formula: where P represents Monte-Carlo simulation iterations and K denotes the number of target signals. (θ p k ,φ p k ) (k = 1, 2 . . . , K) represents the estimate of the k-th angle (θ k , φ k ) for the p-th trial. For this study, P = 500.

Comparison of the DOA Estimation Performance of Different Algorithms
In this subsection, we present the RMSE performance of various algorithms with a CCA, specifically, the PM [42] algorithm, the ESPRIT [43] algorithm, and the proposed TA-MUSIC and E-MUSIC algorithms. It is assumed that there are two DOAs from the directions of (θ 1 , φ 1 ) = (10 • , 20 • ) and (θ 2 , φ 2 ) = (20 • , 30 • ). Figure 7 depicts the RMSE against the SNR (signal-to-noise ratio), where the number of snapshots is L = 200. Figure 8 depicts the RMSE against the number of snapshots, where SNR = 0dB. In Figure 7, it can be seen that estimation performance improved as with increasing SNR, and that the proposed TA-MUSIC algorithm outperforms the others, because it utilizes all the available information, from both the auto-covariance matrix and the mutual-covariance matrix. In addition, Figures 7 and 8 show that the E-MUSIC algorithm, which uses a successive scheme to estimate the signals' direction of arrival, can achieve prominent DOA estimation performance, along with the high computational efficiency demonstrated in Figures 2 and 3.

RMSE with a Varying Number of Sensors
In this subsection, we illustrate how DOA estimation performance changes with the number of snapshots for the proposed E-MUSIC algorithm with CCA geometry which consists of two subarrays. The first subarray, M 1 × N 1 × J 1 , was simulated with various geometries, where the number of sensors was M 1 = N 1 =J 1 = [4, 5,7]. The other subarray remained set to M 2 × N 2 × J 2 = 3 × 3 × 3 sensors throughout. Figure 9 depicts RMSE performance with increasing SNR levels for the various-sized sensor arrays, where L = 200. Figure 10 shows RMSE performance as the number of snapshots increases, where SNR = 0dB. As can be seen in Figure 9, DOA estimation performance improved as the SNR increased, due to the decrease in noise. Figure 10 clearly shows that the DOA estimation performance of the E-MUSIC algorithm improved with the number of snapshots.

RMSE with a Varying Number of Sensors
In this subsection, we illustrate how DOA estimation performa

RMSE with a Varying Number of Sensors
In this subsection, we illustrate how DOA estimation performan

Resolution Performance
If we assume that two target signals are closely located and im array, and define the two signals as 1

Resolution Performance
If we assume that two target signals are closely located and im array, and define the two signals as 1 Figure 10. RMSE of various geometries against the number of snapshots.

Resolution Performance
If we assume that two target signals are closely located and impinge on the CCA array, and define the two signals as (θ 1 ,φ 1 ) and (θ 2 ,φ 2 ), then these two signals are detected if θ −θ < |θ 1 − θ 2 |/2 and φ −φ < |φ 1 − φ 2 |/2. We selected two closely located signals of (θ 1 , θ 2 ) = (10 • , 11 • ) and (φ 1 , φ 2 ) = (30 • , 32 • ) to study the estimation probability of the proposed algorithms. Figures 11 and 12 show the estimation probability of the proposed TA-MUSIC and E-MUSIC algorithms and the conventional ESPRIT and PM algorithms according to SNR and number of snapshots, respectively. As depicted in these two figures, the resolution performance improved with increasing SNR and number of snapshots. Specifically, from Figure 11, it can be seen that the DOA estimation improved with the improved SNR, and Figure 12 clearly shows that the DOA estimation performance improved with the number of snapshots. It can also be seen that the conventional ESPRIT and PM algorithms presented almost identical DOA estimation performances, but the proposed TA-MUSIC and E-MUSIC algorithms demonstrated far better estimation performances.
Sensors 2022, 21, x FOR PEER REVIEW Figure 11. Resolution against SNR for various algorithms.

Conclusions
In this paper, we propose using a TA-MUSIC algorithm for 2D DOA e a CCA configuration. Unlike the conventional MUSIC algorithms which

Conclusions
In this paper, we propose using a TA-MUSIC algorithm for 2D DOA est a CCA configuration. Unlike the conventional MUSIC algorithms which on formation from the auto-covariance matrix, the TA-MUSIC algorithm em

Conclusions
In this paper, we propose using a TA-MUSIC algorithm for 2D DOA estimation with a CCA configuration. Unlike the conventional MUSIC algorithms which only utilize information from the auto-covariance matrix, the TA-MUSIC algorithm employs all the available information from both this and the mutual-covariance matrix. This enables the full achievement of the DOF and improves DOA estimation performance. Moreover, in order to reduce computational complexity, we propose using the computationally efficient E-MUSIC for 2D DOA estimation. This algorithm exploits the ESPRIT algorithm to obtain an initial estimate. The simulation results verified that the proposed algorithms demonstrate superior DOA estimation performance compared to the conventional algorithms.